setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
library(openxlsx)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210316LIBD.eQTL处理/53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
#file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

sci=read.csv("E:/0 公共数据库差异情况/2018 Science supp/aat8127_Table_S1a_DGE.csv",header=T)		###DEG数据合并
sci_scz=sci[sci$SCZ.fdr<0.1,]	#6324 sig gene
sci_bd=sci[sci$BD.fdr<0.1,]		#2010 sig gene
sci_scz=data.frame(symbol=sci_scz$gene_name,SCZ.log2FC=sci_scz$SCZ.log2FC,SCZ.2018sci.fdr=sci_scz$SCZ.fdr)
sci_bd=data.frame(symbol=sci_bd$gene_name,BD.log2FC=sci_bd$BD.log2FC,BD.2018sci.fdr=sci_bd$BD.fdr)
#sci_scz$group.2018sci.scz="SCZ"
#sci_bd$group.2018sci.BD="BD"

cmc.nosva=read.xlsx("E:/0 公共数据库差异情况/CMC/DEG_nosva.xlsx",sheet=1)
cmc.nosva=cmc.nosva[cmc.nosva$adj.P.Val<0.1,]	#1265 sig gene
cmc.nosva=cmc.nosva[!cmc.nosva$MAPPED_genes==".",]
cmc.nosva=data.frame(symbol=cmc.nosva$MAPPED_genes,cmc.logFC=cmc.nosva$logFC,cmc.nosva.fdr=cmc.nosva$adj.P.Val)
cmc.nosva$cmc.nosva.group="CMC.nosva"

#deg=merge(sci_scz,sci_bd,by="symbol",all=T)
#deg=merge(deg,cmc.nosva,by="symbol",all=T)

con=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201207/con_genotype.random3.anno.hg19_multianno.csv",header=T)
conid=unique(unlist(strsplit(con$Gene.refGene,";")))
conid=conid[!conid=="."]

c=length(intersect(as.character(cmc.nosva$symbol),as.character(conid)))
d=length(conid)-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.150K",c,NA,NA)

fileaid=unique(unlist(strsplit(filea$Gene.refGene,";")))
fileaid=fileaid[!fileaid=="."]
a=length(intersect(as.character(fileaid),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(fileaid)-a,d),nrow = 2,byrow = T))
result[2,]=c("117012",a,tmp$estimate,tmp$p.value)

filea1id=unique(unlist(strsplit(filea1$Gene.refGene,";")))
filea1id=filea1id[!filea1id=="."]
a=length(intersect(as.character(filea1id),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(filea1id)-a,d),nrow = 2,byrow = T))
result[3,]=c("61725",a,tmp$estimate,tmp$p.value)

fileid=unique(unlist(strsplit(file$Gene.refGene,";")))
fileid=fileid[!fileid=="."]
a=length(intersect(as.character(fileid),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(fileid)-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

file1id=unique(unlist(strsplit(file1$Gene.refGene,";")))
file1id=file1id[!file1id=="."]
a=length(intersect(as.character(file1id),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(file1id)-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

file2id=unique(unlist(strsplit(file2$Gene.refGene,";")))
file2id=file2id[!file2id=="."]
a=length(intersect(as.character(file2id),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(file2id)-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

file3id=unique(unlist(strsplit(file3$Gene.refGene,";")))
file3id=file3id[!file3id=="."]
a=length(intersect(as.character(file3id),as.character(cmc.nosva$symbol)))
tmp=fisher.test(matrix(c(a,c,length(file3id)-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"./20210317.DEG富集/CMCenrich.con.150K.random3.csv")

library(ggplot2)
scz=read.csv("./20210317.DEG富集/SCZenrich.con.150K.random1.csv",head=T)
bd=read.csv("./20210317.DEG富集/BDenrich.con.150K.random1.csv",head=T)
cmc=read.csv("./20210317.DEG富集/CMCenrich.con.150K.random1.csv",head=T)

library(ggrepel)
library(ggplot2)
scz$group="SCZ"
bd$group="BD"
cmc$group="CMC"
deg=rbind(scz[6,],bd[6,],cmc[6,])

names(deg)=c("X","Term","counts",names(deg)[4:6])
deg$log2Pvalue=-log2(deg$P.value)
ggplot(deg,aes(log2Pvalue,OR))+geom_point(aes(size=3,color=log2Pvalue))+
  scale_color_gradient(low = "#3C5488", high = "#E64B35")+
  theme_light(base_size = 15)+geom_vline(xintercept = -log2(0.05),color="red",linetype="dashed")+
  geom_text_repel(data = deg,aes(label = paste(group,counts,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+
  xlab("-log2(P value)")+ theme(legend.position = "none")
